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^SPECIFI CATION 
~, TITLE 

" METH ?i?,^¥^ B ^ FOR FAST PROCESSING OF MEASUREMENT data 
WITH A PLURALITY OF INDEPENDENT RANDOM SAMPLES 
BACKGROUND nc TH £ INVPMTIOM 

Field of the Invention 

The present invention concerns a method as well as a device to process 
measurement data formed by data sets with a number of independent random 
samples originating from temporally successive measurements, wherein for each 
independent random sample comprised in the data set, a measured time curve (that 
arises from the measurement data of the random sample) is compared using the 
general linear model with the time curve of at least one model function comprised in 
a model matrix, in order to test the incidence of specific patterns in the signal curve. 

More particularly, the present invention concerns the processing of 
measurement data that are acquired from a subject volume with a measurement 
method of functional Imaging, in particular with functional magnetic resonance 
tomography (fMRI), and that are composed of volume data sets originating from a 
number of temporally successive measurements, with, for each volume element 
acquired in the subject volume, a temporal signal curve that arises from 
measurement data of the volume element being compared with the time curve of at 
least one model function. 
Description of the Prior Art 

In particular in the field of medical technology and medical research, a need 
exists to acquire information about the brain activity in human and animal organs. 
The neural activity causes an increase of the blood flow in active brain areas, 
resulting in a decrease in the blood-deoxyhemoglobin concentration. 
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Deoxyhemoglobin is a paramagnetic material that reduces the magnetic field 
homogeneity and therefore can be shown with the aid of magnetic resonance 
techniques, since it accelerates the T 2 * signal relaxation. 

Localization of brain activity is enabled by the use of a functional imaging 
method that measures the change of the MR signal relaxation with a time delay. The 
biological effective mechanism is known in the literature by the name BOLD effect 
(Blood Oxygen Level Dependent effect), 

A fast magnetic resonance imaging enables the examination of the bold effect 
in vivo, dependent on activation states of the brain. In functional magnetic 
resonance tomography, magnetic resonance exposures of the subject volume to be 
examined, for example of the brain of a patient, are made in short temporal intervals. 
By comparison of the signal curve (measured by functional imaging) for each volume 
element of the subject volume with the time curve of a model function, a stimulus- 
specific neural activation can be detected and spatially localized. A stimulus can be, 
for example, a somatic sensory, acoustic, visual or olfactory stimulus, as well as a 
mental or motor task. The model function or the model time series specifies the 
expected signal change of the magnetic resonance signal in the course of neural 
activation. By using faster magnetic resonance techniques, such as, for example, 
the echo-planar method, smaller temporal intervals between the individual 

» 

measurements can be realized. 

In many multivariate statistical analyses, a model known as the general linear 
model (GLM) is used for the comparison of the measured signal curve with the time 
curve of a model function. With the general linear model, It Is determined which 
Ifnear combination of the model functions best approximates the measurement data 
series. Furthermore, it can be calculated for each model function how significantly 
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the measurement data of the null hypothesis of no contribution of the respective 
model function contradicts the measurement data series. The general linear model 
is used in many fields, such as, for example, physics and sociology, to analyze 
measurement data. It is in particular also suitable to analyze time series as they are 
measured in functional magnetic resonance imaging (fMRI). By using the general 
linear model, it can be analyzed whether the measured time series show a pattern 
that corresponds to local neural activity. In addition to this pattern, however, the time 
series frequently also show other characteristics (such as, for example, drifts or other 
effects) that can likewise be modeled in the framework of the general linear model. 
This enables a better analysis of the measurement data than, for example, a West or 
correlation method, Thus, for example, it is also possible with the general linear 
model to analyze in parallel a number of effects in the brain. Group statistics about a 
number of test persons are also possible. Further application possibilities of the 
general linear model are found, for example, in "Human Brain Function" by R. 
Frackowiak et al M Academic Press, 

In the processing of measurement data that are acquired from a subject 
volume with the method of functional magnetic resonance tomography, it has been 
necessary until now to load in the main storage of a computer the overall 
measurement data that comprises a number of volume data sets originating via 
temporally successive measurements. Subsequently, for each volume element of 
the measured subject volume, the signal curve or, respectively, the time series must 
be extracted from these measurement data and be compared with the model 
function. In the known Implementation of the general linear model in the freely 
available SPM software (Wellcome Department of Cognitive Neurology; University of 
London; published under Gnu Public License; http://www.fil.ion.ucl.ac.uk/spm/), it is 
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likewise necessary to load into the main storage of the computer the complete data 
set to be analyzed. The complete data set, in long fMRI studies - possibly also 
spanning a number of test persons - can include several hundred megabytes, up to 
gigabytes, of data. The values to be analyzed that belong to a time series of 
measurement data are extracted, and the general linear model is directly calculated. 

This conventional data processing therefore leads to a significant main 
storage requirement Since the measurement data typically exist in the storage of 
the computer volume by volume, corresponding to the successive measurements as 
a number of volume data sets, the use of this known technique also leads to very 
long computing times, since the individual time series must be collected together 
over very large ranges of the loaded data. 

Presorting of the data could reduce this computing time, however this requires 
in turn a considerable computing time and additional storage requirements. 
Moreover, a sorting event can not be finished until the end of the measurement, 
since the time series only then exists in full. The actual calculations thus can be 
started only at the end of the measurement. 

SUMMARY OF THE INVENTION 

An object of the invention is to provide a method and a device for fast and 
storage-saving processing of measurement data that are formed of a number of data 
sets with a plurality of independent random samples, originating from temporally 
successive measurements. The method should in particular enable fast and 
storage-saving calculation of the general linear model for the case of a number of 
independent random samples, for example measurement data series as they occur 
in fMRI. 
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This object is achieved in accordance with the invention for the processing of 
measurement data that is composed of a number of data sets with a number of 
independent random samples, originating from temporally successive 
measurements. 

In the following, the method and its advantages are explained (without 
limitation to other applications) using the processing of measurement data that are 
acquired from a subject volume with a measurement method of functional imaging, in 
particular with the method of functional magnetic resonance tomography, and that 
are formed of a number of volume data sets originating from temporally successive 
measurements. In this method, for each volume element acquired in the subject 
volume, a temporal signal curve is determined from the measurement data of this 
volume element, and at least one model function obtained in a model matrix G using 
the general linear model is compared with the time curve in order to test the 
incidence of specific characteristics or patterns in the signal curve. The present 
method is primarily distinguished in that calculations necessary for the comparison 
from the measurement data in a sequence of the of the data sets, or volume data 
sets originating from the time sequence of the measurements, are implemented for 
all relevant measurement data of a volume data set and are stored as intermediate 
results. The intermediate results from the directly preceding volume data set are 
updated with the new calculations, such that at any time intermediate results can be 
calculated, and after the one-time cycle of all volume data sets an end result exists 
that allows an assertion about the incidence of the characteristics or patterns in the 
signal curve. The square of an error vector necessary for the end result is Obtained 
from the difference of the square of a measurement value vector formed from the 
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measurement data and the square of a model vector calculated from th© model 
function- 

In contrast to the known method of the prior art, in the present method the 
measurement data are not cycled for the comparison according to time series or 
volume elements, but rather volume-by-volume. After implementation of the 
calculations with the measurement data and storage of an intermediate result, each 
volume data set can again be discarded from the main storage of the data 
processing system. For this reason, it is no longer necessary to load the entirety of 
the measurement data, meaning all volume data sets of the measurements, into the 
working storage of the data processing system. Rather P it is sufficient merely to load 
the measurement data of an individual volume data set For this reason, the present 
method can be implemented efficiently and quickly, saving much storage space, 
since respectively only small data ranges must be accessed in the calculation. The 
volume data sets must hereby overall only be cycled once, volume-by-volume. 

In the solution of the above object, consideration has been made that the 
present computer architecture is very slow in the calculation of large data ranges, in 
particular in the gigabyte range, however it can implement small data ranges with 
noticeably higher speed, With the present method, in which only small data ranges 
are always used for calculation volume-by-volume, a noticeably higher calculation 
speed thus can be achieved. Furthermore, the storage requirement that is 
necessary for the processing of the data is no longer dependent in the present 
method on the number of the volume data sets (meaning the number of the 
measurements), but rather is substantially dependent only on the size of the 
measured subject volume. A further advantage of the present method is that all 
calculations can already be implemented during the measurement. Real-time 
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applications as well as the display of an always-current intermediate result are also 
thereby possible. Furthermore, the model function in the present method can first be 
defined during the measurement, and does not have to be known a priori. The 
method can be implemented in parallel very efficiently, such that a number of 
computers or processors can be used to accelerate the calculations. 

The present method is explained again as an example using an fMRI data set, 
however it can naturally also be applied in the same manner to data sets of other 
functional imaging methods or other measurement data with independent random 
samples, corresponding to the measurement data series of the individual volume 
elements in the fMRI, 

In functional magnetic resonance tomography, the subject volume, for 
example the brain of a test person, is repeatedly sampled three-dimensionally at 
small temporal intervals, and from the raw data the desired image information is 
reconstructed in voxels, meaning volume elements, via a Fourier transformation. A 
corresponding measurement value is hereby associated with each volume element. 
The overall measurement is composed of a number of volume data sets that form 
the measurement data of the measurements temporally following in quick 
succession. Each volume data set is composed of the measurement data, acquired 
at a specific point in time, of the individual voxel of the subject volume. For the 
recognition of activation states of the measured brain, the signal curves resulting 
from the temporally successive measurement values of each voxel are compared in 
order to be able to establish the degree of a coincidence. Due to the use the 
widespread general linear model (GLM) according to the present method, this 
comparison can be implemented with a great saving of storage space. In the 
general linear model, a design matrix or model matrix G is generated with one or 
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more model functions from which, in connection with the measurement data of each 
voiume element, values (in particular a t-statistic) can be caicuiated, from which the 
degree of coincidence with the model functions of the model matrix is visible. 

In the implementation of the method using the general linear model, the 
volume data sets of the measurements are only cycled once. After cycling each data 
set with the thereby-implemented calculations, an intermediate result is stored that is 
updated upon cycling of the directly ensuing volume data set via the new 
calculations - due to the new measurement data - and is stored as an intermediate 
result of this new volume data set. By this volume-by-volume updating of the 
intermediate results, it is possible to obtain a result at any time from which, for 
example, the t-statistic can be calculated. The square of the error vector required for 
this is obtained from the difference of the square of the measurement value vector 
formed from the measurement data and the square of the model vector calculated 
from the model function. This simple calculation is possible since the error vector 
and the model vector are perpendicular to one another, due to the least squares 
condition, and when added yield the measurement value vector. For the calculation 
of the square of the error vector, therefore only the length of the measurement value 
vector and the length of the model vector must be determined. The square of the 
error vector can then be calculated by the set of Pythagoras. The model vector 
thereby comprises the values calculated from the model function that should 
approach as best possible the measurement values. 

The present method is thus based on an update event in which the volume 
data sets are only cycled once. At the end of the update event, the data necessary 
for the calculation of the t-statistic from the measurement values are all available. 
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in one embodiment, the update event, the cycle throu S h the volume date sets 
can ensue simultaneous* with the measurement, with the respective ca.cula.ion 
being implemented directly after obtaining a new volume data set. In this 
embodiment, the update event is therefore already finished with the end of the data 
acquisition. 

Given sufficient computing speed, a real-time application of the present 
method ,s also possible, in that during the measurement regular intermediate 
evaluations are undertaken in which me update event is interrupted and the 
calculation of the t-statistic is implemented. The update event can subsequently be 
continued for the volume data sets measured in the following time period. 

Ih the present method, buffers are necessary that, for the use of the general 
linear model, must exhibit a size that satisfy the following equation: 
Size of the storage = number of the volume elements x (2 x number of the model 
functions + 2) x 4 bytes. 

The four bytes as a storage space for a floating-point number of the float data 
type have turned out to be sufficient for the implementation of the present method. 
The number of the volume elements can, for example, comprise 64 x 64 x 32 values. 
This corresponds to the size of a typical sample data set of a subject volume in 
functional magnetic resonance tomography. 

The present method can be used, for example, in a computer of a magnetic 
resonance scanner, such that a fast postprocessing of fMRI data is possible. Via 
suitable implementation of the method, the representation of the comparison results 
can a , so ensue in real-time, as ,ong as the data sets do not comprise too many 
volume elements and the number of the model functions is not chosen too high. In 
particular, the calculation in reaUime is not limited by the number of the volume data 
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sets, since the calculation time is constant, independent of tho number of , ha already 
added data sets. 

DESCRIPT ION OF THE DRAWINGS 

The figure is a flow chart of an exemplary embodiment of the inventive 
method. 

DESCRIPTION OF TUP P REFERRED EMBODIMENTS 

In the inventive method, first the model matrix G is defined, that is composed 
of a number of columns with model functions. The time series in the present 
invention merely contains twenty measurements. For each volume element, twenty 
measurement values are thereby obtained, whereby m = 0 ... rna*1 and max = 
20. As a model matrix, a matrix wim five model functions is selected, meaning a 
matrix with five columns and twenty rows. 

In the first column is the model function of the neural activation that for 
example can be formed by: 

G = (11 1110001 1001 1001 100) T 
In the second column, a constant ration is represented in the measurement data: 

Gm,i ~ 1. 

The further three columns serve as high-pass filters, meaning they model 
slow drifts during the measurement on a cosine basis [cosinus basis]. These model 
functions can. for example, be selected as 

(=) fiS — (^] 

The model matrix G„j thus comprises the collection of model functions that 
should be adapted to the measurement data Xm. 
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In the implementation of the comparison m h, 
model function, , 00n V-*»< of the measurement date with these 

vei rr: quares e ~ ^ ~ » A 
si 2 ' ; : fram whioh ,he — - 

calculated in the following manner: 



B=(G T G)-i G T x. 



vector b in „. presen , examp(e fe 
-spond.o to the „ „ „ ^ ^ _ ^ ^ ■ 

«. -ue represent ho w the co^pondin, ^ ^ 
order to approximate the measurement data. 

aesenloTr" 7 * — <° — - 

zrrr. qua,i,y of ,he es8ma,e - fw •* - va,ue ■ is - 

results in the following manner: 



Vmax-rang(G) 



whereby e - G ' b - x. 



-nall, what is Wr, as fhe ,-sta,i ste „ catenated by . for each volume 
*™nt. ootalnln, a value , havin 9 a ma g n„ude that is . measure fo r the ee 9 ree of 

Z7 r of ,he ~ em date series ~ - - «•— 

, F 7 here as an examp,e ' * is *~ to — - — 
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specifies the model function of interest. In the present case, this contrast c can be 
given by the following value: c = ^ 0000) T . 

The t-value then results in the following manner: 



t ■= 0 ' k 



In the implementations implemented until now in the prior art, the time series 
of each voxel is extracted in succession from the fMRI data set and undergoes the 
analysis above. However, for this it is necessary for the calculation to 
simultaneously load all measurement data of all volume data sets into the working 
storage of a computer. In addition to the high storage space requirement, this also 
leads to a correspondingly slow calculation speed. 

In contrast to this, with the present invention as shown in the figure in an 
embodiment, both the storage space requirement can be reduced and the 
calculation speed can be increased. In the inventive method, in the update event the 
fMRI data set is cycled volume by volume, meaning one volume data set after the 
other, and the respective calculations are implemented in each volume data set for 
each volume element. 

This can ensue in the following manner. In the update event, in which all 
volume data sets are cycled once, the vector GTx is updated for each new volume 
data set and each voxel, whereby respectively only the new line of the model matrix 
and the new measurement value are accessed: 

GTx, = GTXi + G m> , ■ x^ 



12 



NOV. 14. 2003 11 :54PM NO. 188 P. 15/28 

GTX hereby designates the vector G T ■ x. This calculation Is implemented in 
each new volume data set for each volume element, and can also be distributed 
voxel by voxel onto various computers or, respectively, processors. The value 
obtained from the preceding volume data set for the Individual components of the 
vector GTx is hereby updated. Of course, this value must be set to zero in the 
implementation of the calculation for the very first volume data set. 

The implementation of this calculation for each volume data set implies that a 
sufficient buffer must be present for the respective interim value of GTx that 
corresponding to the number of the volume elements in the volume data set 
multiplied by the number of the model functions in the model matrix. 

Furthermore, in this cycle of the data sets for each new volume data set and 
each voxel, the square XX of the measurement value vector us updated, with only 
the new measurement value is accessed: 

XX = XX + (x m ) 2 

At the same time, in this cycle the matrix GTG (corresponds to G T ■ G) can 
also be update by following instruction: 

GTG ilk = GTG,, k + G m ,, ■ G m>k . 

Naturally, this calculation (since it involves no measurement data) also can be 
implemented in a step before the cycle, as indicated dashed in the figure, or after 
this cycle. However, in a model matrix whose model functions are first defined or 
changed in the course of the measurement, or when intermediate results should be 
calculated, the implementation is necessary according to the above instruction, 
meaning the volume-by-volume calculation during the cycle. 
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After obtaining the matrix GTG, this is inverted. This can ensue, for example, 
by means of an LU decomposition, since this matrix is a small, symmetric, real 
matrix, in the present example a 5 x 5 matrix. Naturally, however, other inversion 
algorithms can also be used. Finally, the vector b is calculated for each volume 
element, which results by 

b = GTG' 1 • GTx 

After implementation of the update event, the values for the vector b and the 
pseudo-inverse (G T ' G)" 1 are thus established. 

To calculate the value for the square of the error vector e (designated in the 
following as EE), used is made of the fact that the part of the measurement value 
vector {designated as model vector M) representable via the model, and the vector 
E, are perpendicular to one another, due to the least squares condition, and when 
added result in the measurement value vector X. Thus only the length of the 
measurement value vector and the length of the model vector must be determined in 
order to be able to calculate from this the set of Pythagoras EE. 
Thus: 

M - G • b Model matrix = M ■ M 

F 

MM*(G b) (G b) = q q* £q, .q, 

with: q, -^O^-bj 

The use of which results in: 
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max -J sp-l sp-l 
i-0 >0 k=0 



i=0 



ZI(G,. i 'G l ,)'b J -b, 



The inner parentheses represent the matrix GTG" 1 that was already 
calculated. Thus: 

MM = MM + GTGj.k 1 bj • b k 
such that MM can be calculated quickly, and leads directly to EE: 

EE = XX - MM. 

For the previously implemented calculations: 

i = 0 sp-1 
k = 0 sp-1 
j =» 0 .. sp-1 
m = 0 .. max-1 

whereby max corresponds to the number of the measurements or, respectively, 
volume data sets and sp corresponds to the columns of the model matrix G. 

After each addition of a volume data set, or also at the end of the 
measurement, the t-statistic cannot be calculated for this comparison. In the present 
example, a scaling factor Scale is hereby first calculated according to the following 
instruction: 

Scale = Scale + contrast increase ■ GTGjj ■ q , 



whereby c corresponds to the contrast vector. 
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This calculation is implemented for each volume element. Furthermore, the 
a 2 is calculated for each volume element, which results according to the following 
known instruction: 



a 2 

max- rang(G) 

From this value, what is known as the t-value can be separately calculated for 
each voxel by 

t + Ci ■ b| . 

This t is finally divided again by Eq, in order to obtain for each voxel the final t-value 
that, for example, can be visualized by superimposition on a typical magnetic 
resonance exposure of the subject volume. The final calculations can likewise by 
distributed voxel by voxel to different computers or. respectively, processors 
calculating in parallel, in order to accelerate the calculation. 

Furthermore, the figure shows the possibility, indicated by the dashed arrows, 

to calculate and to output an intermediate result during the cycle of the data sets. 

Although modifications and changes may be suggested by those skilled In the 

art, it is the intention of the inventor to embody within the patent warranted hereon all 

changes and modifications as reasonably and properly come within the scope of his 

contribution to the art. 
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